#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''
@Project ：V2 
@File    ：DB1.py
@IDE     ：PyCharm 
@Author  ：郭星
@Date    ：2025/9/28 18:24 
'''
import numpy as np
import random
from scipy.spatial.distance import cdist
from collections import defaultdict


# -------------------------- 1. 问题参数配置 --------------------------
class ProblemConfig:
    def __init__(self):
        # 流水线与缓冲区数量
        self.n_lines = 3  # 3条流水线（L1-L3）
        self.m_buffers = 8  # 每条流水线8个缓冲区（B1-B8）
        # 缓冲区约束
        self.C_max_total = 1500  # 总容量上限
        self.C_min_single = 1  # 单个缓冲区最小容量
        self.C_max_single = 200  # 单个缓冲区最大容量
        # 传输成本（-表示不可传输）
        self.transfer_cost = np.array([
            [np.inf, np.inf, np.inf],  # B1
            [5, 5, np.inf],  # B2
            [np.inf, 4, 4],  # B3
            [5, 5, np.inf],  # B4
            [np.inf, 4, 4],  # B5
            [5, 5, np.inf],  # B6
            [6, np.inf, 6],  # B7
            [np.inf, np.inf, np.inf]  # B8
        ])
        # 产品收益
        self.product_profit = [4, 4.5, 4, 5, 3.5]
        # 算法参数
        self.pop_size = 60  # 种群规模
        self.max_incnum = 5  # 滞化指数阈值
        self.Er = 0.15  # 子群迁移比例
        self.attr0 = 1.0  # 初始吸引度
        self.gamma = 0.5  # 亮度衰减系数
        self.T0 = 99  # 模拟退火初始温度
        self.T_decay = 0.99  # 温度衰减率
        self.K_max = 0.15  # 高斯扰动上限
        self.K_min = 0.05  # 高斯扰动下限


# -------------------------- 2. 核心工具函数 --------------------------
def non_dominated_sort(population, obj_values):
    """Pareto非支配排序：返回非支配解集与支配等级"""
    n = len(population)
    if n == 0:  # 处理空种群情况
        return [], [], []

    dominated_count = [0] * n  # 被多少个体支配
    dominates_list = [[] for _ in range(n)]  # 支配的个体列表
    fronts = [[]]  # 支配等级（front[0]为最优前沿）

    # 计算支配关系
    for i in range(n):
        for j in range(n):
            if i == j:
                continue
            # 个体i支配个体j
            if (obj_values[i][0] >= obj_values[j][0] and obj_values[i][1] >= obj_values[j][1]) and \
                    (obj_values[i][0] > obj_values[j][0] or obj_values[i][1] > obj_values[j][1]):
                dominates_list[i].append(j)
            elif (obj_values[j][0] >= obj_values[i][0] and obj_values[j][1] >= obj_values[i][1]) and \
                    (obj_values[j][0] > obj_values[i][0] or obj_values[j][1] > obj_values[i][1]):
                dominated_count[i] += 1

        # 初始前沿（无支配个体）
        if dominated_count[i] == 0:
            fronts[0].append(i)

    # 生成后续前沿
    current_front = 0
    while current_front < len(fronts) and len(fronts[current_front]) > 0:
        next_front = []
        for i in fronts[current_front]:
            for j in dominates_list[i]:
                dominated_count[j] -= 1
                if dominated_count[j] == 0:
                    next_front.append(j)
        current_front += 1
        if next_front:
            fronts.append(next_front)

    # 提取非支配解集（最优前沿）
    archive_indices = fronts[0] if fronts else []
    archive = [population[i] for i in archive_indices]
    archive_objs = [obj_values[i] for i in archive_indices]
    return archive, archive_objs, fronts


def crowding_distance(obj_values, front_indices):
    """计算拥挤度：修复索引越界问题，使用前沿内部相对索引"""
    n = len(front_indices)
    if n <= 1:  # 空前沿或只有一个个体时，直接返回空距离或无穷大
        return [np.inf] * n if n == 1 else []

    distance = [0.0] * n  # 长度为当前前沿的个体数量
    obj_dim = len(obj_values[0])

    for dim in range(obj_dim):
        # 关键修复：对前沿内部的索引进行排序（0到n-1），而不是种群全局索引
        # 按当前目标维度排序前沿内的个体
        sorted_indices = sorted(range(n), key=lambda i: obj_values[front_indices[i]][dim])

        # 边界个体拥挤度设为无穷大（使用前沿内部索引）
        distance[sorted_indices[0]] = np.inf
        distance[sorted_indices[-1]] = np.inf

        # 计算中间个体的拥挤度
        # 获取当前维度的最小值和最大值（用于归一化）
        min_val = obj_values[front_indices[sorted_indices[0]]][dim]
        max_val = obj_values[front_indices[sorted_indices[-1]]][dim]

        if max_val - min_val == 0:
            continue  # 目标维度无差异，跳过

        # 计算中间个体的拥挤度
        for i in range(1, n - 1):
            # 前一个个体的目标值
            prev_val = obj_values[front_indices[sorted_indices[i - 1]]][dim]
            # 后一个个体的目标值
            next_val = obj_values[front_indices[sorted_indices[i + 1]]][dim]
            # 累加拥挤度（使用前沿内部索引）
            distance[sorted_indices[i]] += (next_val - prev_val) / (max_val - min_val)

    return distance


# -------------------------- 3. C-MOHFA算法实现 --------------------------
class C_MOHFA:
    def __init__(self, config):
        self.config = config
        self.population = []  # 种群：(容量向量C, 传输率向量TS)
        self.obj_values = []  # 目标值：(总收益W, 平均吞吐率theta)
        self.archive = []  # 非支配解集
        self.archive_objs = []  # 非支配解的目标值
        self.buffer_archive = defaultdict(list)  # 缓冲区档案集
        self.incnum = 0  # 滞化指数
        self.iter = 0  # 迭代次数

    def initialize_population(self):
        """多策略初始化种群：随机+启发式+生产率"""
        pop_size = self.config.pop_size
        n, m = self.config.n_lines, self.config.m_buffers
        C_max_total = self.config.C_max_total

        # 1. 随机初始化（30%）
        rand_num = int(pop_size * 0.3)
        for _ in range(rand_num):
            # 随机生成容量（满足总容量约束）
            C = np.random.randint(self.config.C_min_single, self.config.C_max_single + 1, (n, m))
            C = C * (C_max_total // C.sum())  # 缩放到总容量约束内
            if C.sum() < C_max_total:
                C[np.random.randint(n)][np.random.randint(m)] += C_max_total - C.sum()
            # 随机生成传输率（[0,1]）
            TS = np.random.uniform(0, 1, (n, m))
            # 不可传输的缓冲区设为0
            for j in range(m):
                for i in range(n):
                    if self.config.transfer_cost[j][i] == np.inf:
                        TS[i][j] = 0.0
            self.population.append((C, TS))

        # 2. 启发式初始化（30%）：高斯分布权重分配
        heu_num = int(pop_size * 0.3)
        alpha = [0.52, 0.49, 0.47]
        beta = [0.21, 0.19, 0.23]
        for _ in range(heu_num):
            C = np.zeros((n, m), dtype=int)
            for i in range(n):
                mu = alpha[i] * m
                sigma = beta[i] * m
                weights = np.random.normal(mu, sigma, m)
                weights = np.maximum(weights, 0.1)  # 避免权重为0
                weights = weights / weights.sum()  # 归一化
                # 分配容量
                line_C = np.round(weights * (C_max_total / n)).astype(int)
                if line_C.sum() < (C_max_total / n):
                    line_C[np.argmax(weights)] += (C_max_total / n) - line_C.sum()
                C[i] = line_C
            # 传输率：基于权重生成
            TS = np.zeros((n, m))
            for i in range(n):
                weights = C[i] / C[i].sum()
                for j in range(m):
                    if self.config.transfer_cost[j][i] != np.inf:
                        TS[i][j] = weights[j] * np.random.uniform(0.5, 1.0)
            self.population.append((C, TS))

        # 3. 生产率初始化（40%）：基于设备产能分配
        prod_num = pop_size - rand_num - heu_num
        capacity_mean = np.array([
            [17.5, 10.5, 19.5, 11.0, 17.0, 19.0, 13.5, 11.0, 16.0],  # L1
            [16.5, 18.5, 10.5, 19.0, 10.5, 13.5, 20.0, 15.0, 13.5],  # L2
            [19.5, 14.0, 16.5, 12.5, 18.5, 12.0, 10.0, 16.0, 15.0]  # L3
        ])
        buffer_demand = np.abs(capacity_mean[:, :-1] - capacity_mean[:, 1:])
        for _ in range(prod_num):
            demand_total = buffer_demand.sum()
            C = np.round(buffer_demand * (C_max_total / demand_total)).astype(int)
            if C.sum() < C_max_total:
                C[np.unravel_index(np.argmax(buffer_demand), buffer_demand.shape)] += C_max_total - C.sum()
            # 传输率：需求高的缓冲区传输率高
            TS = np.zeros((n, m))
            for i in range(n):
                demand_norm = buffer_demand[i] / buffer_demand[i].sum()
                for j in range(m):
                    if self.config.transfer_cost[j][i] != np.inf:
                        TS[i][j] = demand_norm[j] * np.random.uniform(0.5, 1.0)
            self.population.append((C, TS))

        # 初始化目标值
        self.obj_values = [self.evaluate_objective(C, TS) for C, TS in self.population]
        # 初始非支配解集
        self.archive, self.archive_objs, _ = non_dominated_sort(self.population, self.obj_values)
        # 更新缓冲区档案集
        for idx, (C, TS) in enumerate(self.archive):
            C_key = tuple(C.flatten())
            self.buffer_archive[C_key].append((TS, self.archive_objs[idx][0]))

    def evaluate_objective(self, C, TS):
        """评估目标函数：修复空切片均值问题"""
        n, m = self.config.n_lines, self.config.m_buffers
        equipment_availability = np.random.uniform(0.85, 0.95, (n, 9))
        line_production = []
        for i in range(n):
            # 设备产能（分流水线）
            if i == 0:
                capacity = np.array([17.5, 10.5, 19.5, 11.0, 17.0, 19.0, 13.5, 11.0, 16.0])
            elif i == 1:
                capacity = np.array([16.5, 18.5, 10.5, 19.0, 10.0, 13.5, 20.0, 15.0, 13.5])
            else:
                capacity = np.array([19.5, 14.0, 16.5, 12.5, 18.5, 12.0, 10.0, 16.0, 15.0])

            effective_capacity = capacity * equipment_availability[i]
            production = np.min(effective_capacity)  # 瓶颈设备决定产量

            # 修复空切片均值问题：当TS[i]全为0时，用0作为默认值
            valid_ts = TS[i][TS[i] > 0]
            share_boost = 0.05 * np.mean(valid_ts) if len(valid_ts) > 0 else 0.0
            production *= (1 + share_boost)
            line_production.append(production)

        # 计算总收益
        product_profit = np.random.choice(self.config.product_profit, n)
        output_profit = sum(line_production[i] * product_profit[i] for i in range(n))
        # 传输成本
        transfer_volume = C * TS
        transfer_cost = sum(
            transfer_volume[i][j] * self.config.transfer_cost[j][i]
            for i in range(n) for j in range(m)
            if self.config.transfer_cost[j][i] != np.inf
        )
        total_profit = output_profit - transfer_cost

        # 计算平均吞吐率
        raw_material = [p * 1.05 for p in line_production]
        line_throughput = [line_production[i] / raw_material[i] for i in range(n)]
        avg_throughput = np.mean(line_throughput)

        return (total_profit, avg_throughput)

    def search_transfer_rate(self, C):
        """传输率搜索：结合档案集+模拟退火+邻域搜索"""
        n, m = self.config.n_lines, self.config.m_buffers
        C_key = tuple(C.flatten())
        T = self.config.T0

        # 优先从档案集获取初始TS
        if C_key in self.buffer_archive:
            ts_list = self.buffer_archive[C_key]
            initial_TS = max(ts_list, key=lambda x: x[1])[0]
        else:
            initial_TS = np.zeros((n, m))
            for i in range(n):
                for j in range(m):
                    if self.config.transfer_cost[j][i] != np.inf:
                        initial_TS[i][j] = np.random.uniform(0.2, 0.8)

        # 模拟退火+邻域搜索优化TS
        current_TS = initial_TS.copy()
        current_obj = self.evaluate_objective(C, current_TS)[0]
        best_TS = current_TS.copy()
        best_obj = current_obj

        while T > 0.1:
            new_TS = current_TS.copy()
            i = np.random.randint(n)
            j = np.random.randint(m)
            if self.config.transfer_cost[j][i] != np.inf:
                new_TS[i][j] = np.clip(new_TS[i][j] + np.random.normal(0, 0.1), 0.0, 1.0)
            else:
                new_TS[i][j] = 0.0

            new_obj = self.evaluate_objective(C, new_TS)[0]
            delta = new_obj - current_obj

            if delta > 0 or np.exp(delta / T) > np.random.random():
                current_TS = new_TS
                current_obj = new_obj
                if current_obj > best_obj:
                    best_TS = current_TS
                    best_obj = current_obj

            T *= self.config.T_decay

        self.buffer_archive[C_key].append((best_TS, best_obj))
        return best_TS

    def firefly_migration(self):
        """萤火虫迁移：修复空前沿导致的索引错误"""
        n_pop = len(self.population)
        n, m = self.config.n_lines, self.config.m_buffers
        new_population = []

        # 计算拥挤度（跳过空前沿）
        _, _, fronts = non_dominated_sort(self.population, self.obj_values)
        crowding_dist = [0.0] * n_pop

        for front in fronts:
            if not front:  # 跳过空前沿
                continue
            dist = crowding_distance(self.obj_values, front)
            for i, idx in enumerate(front):
                crowding_dist[idx] = dist[i]

        # 遍历每个萤火虫（个体）
        for x in range(n_pop):
            Cx, TSx = self.population[x]
            # 选择吸引源y：优先非支配前沿，同前沿选拥挤度低的
            candidates = [y for y in range(n_pop) if y != x]
            if not candidates:  # 避免种群只有一个个体时出错
                new_population.append((Cx, TSx))
                continue

            # 按支配等级排序
            front_level = {}
            for y in candidates:
                for l, f in enumerate(fronts):
                    if y in f:
                        front_level[y] = l
                        break
            candidates.sort(key=lambda y: (front_level[y], crowding_dist[y]))
            y = candidates[0]
            Cy, TSy = self.population[y]

            # 计算欧氏距离
            dist_C = cdist(Cx.flatten().reshape(1, -1), Cy.flatten().reshape(1, -1))[0][0]
            dist_TS = cdist(TSx.flatten().reshape(1, -1), TSy.flatten().reshape(1, -1))[0][0]

            # 计算吸引度
            attr = self.config.attr0 * np.exp(-self.config.gamma * (dist_C + dist_TS))

            # 高斯扰动
            K = self.config.K_max - (self.config.K_max - self.config.K_min) * (self.iter / 100)
            noise_C = np.random.normal(0, K, Cx.shape)
            noise_TS = np.random.normal(0, K, TSx.shape)

            # 更新位置
            new_C = np.round(Cx + attr * (Cy - Cx) + noise_C).astype(int)
            new_TS = np.clip(TSx + attr * (TSy - TSx) + noise_TS, 0.0, 1.0)

            # 修复容量约束
            new_C = np.maximum(new_C, self.config.C_min_single)
            new_C = np.minimum(new_C, self.config.C_max_single)
            if new_C.sum() > self.config.C_max_total:
                excess = new_C.sum() - self.config.C_max_total
                while excess > 0:
                    i = np.random.randint(n)
                    j = np.random.randint(m)
                    if new_C[i][j] > self.config.C_min_single:
                        new_C[i][j] -= 1
                        excess -= 1
            # 修复传输率约束
            for j in range(m):
                for i in range(n):
                    if self.config.transfer_cost[j][i] == np.inf:
                        new_TS[i][j] = 0.0

            new_population.append((new_C, new_TS))

        # 替换旧种群并重新评估
        self.population = new_population
        self.obj_values = [self.evaluate_objective(C, TS) for C, TS in self.population]

    def population_exchange(self):
        """种群信息交换：子群迁移"""
        n_pop = len(self.population)
        if n_pop <= 1:  # 避免种群过小时出错
            return
        exchange_num = max(1, int(n_pop * self.config.Er))  # 确保至少迁移1个个体
        exchange_indices = random.sample(range(n_pop), exchange_num)
        exchange_individuals = [self.population[i] for i in exchange_indices]
        # 移除迁移个体
        self.population = [ind for i, ind in enumerate(self.population) if i not in exchange_indices]
        # 随机分配到种群
        for ind in exchange_individuals:
            pos = random.randint(0, len(self.population) - 1)
            self.population.insert(pos, ind)
        # 重新评估
        self.obj_values = [self.evaluate_objective(C, TS) for C, TS in self.population]

    def run(self, max_iter=100):
        """算法主运行流程"""
        self.initialize_population()
        if not self.archive_objs:  # 处理空解集情况
            print("初始种群为空或无效")
            return None, None, None

        print(
            f"初始非支配解集规模：{len(self.archive)}，最优收益：{max(self.archive_objs, key=lambda x: x[0])[0]:.2f}，最优吞吐率：{max(self.archive_objs, key=lambda x: x[1])[1]:.4f}")

        # 迭代优化
        while self.iter < max_iter and self.incnum < self.config.max_incnum:
            self.iter += 1
            # 1. 传输率重新优化
            new_pop = []
            for C, TS in self.population:
                new_TS = self.search_transfer_rate(C)
                new_pop.append((C, new_TS))
            self.population = new_pop
            self.obj_values = [self.evaluate_objective(C, TS) for C, TS in self.population]

            # 2. 种群信息交换
            self.population_exchange()

            # 3. 萤火虫迁移
            self.firefly_migration()

            # 4. 更新非支配解集
            new_archive, new_archive_objs, _ = non_dominated_sort(self.population, self.obj_values)
            combined_archive = self.archive + new_archive
            combined_objs = self.archive_objs + new_archive_objs
            combined_archive, combined_objs, _ = non_dominated_sort(combined_archive, combined_objs)

            # 判断是否更新
            if len(combined_archive) > len(self.archive) or \
                    (len(combined_archive) == len(self.archive) and max(combined_objs, key=lambda x: x[0])[0] >
                     max(self.archive_objs, key=lambda x: x[0])[0]):
                self.archive = combined_archive
                self.archive_objs = combined_objs
                self.incnum = 0
            else:
                self.incnum += 1

            # 打印迭代信息
            if self.iter % 10 == 0:
                best_profit = max(self.archive_objs, key=lambda x: x[0])[0] if self.archive_objs else 0
                best_theta = max(self.archive_objs, key=lambda x: x[1])[1] if self.archive_objs else 0
                print(
                    f"迭代{self.iter:2d}：非支配解集规模={len(self.archive)}，最优收益={best_profit:.2f}，最优吞吐率={best_theta:.4f}，滞化指数={self.incnum}")

        # 输出最终结果
        print("\n算法终止！")
        if not self.archive_objs:
            print("未找到有效解")
            return None, None, None

        final_best = max(self.archive_objs, key=lambda x: x[0] * x[1])
        final_C, final_TS = self.archive[self.archive_objs.index(final_best)]
        print(f"最终最优解：总收益={final_best[0]:.2f}，平均吞吐率={final_best[1]:.4f}")
        print(f"最优缓冲区容量配置（3条线×8个缓冲区）：\n{final_C}")
        print(f"最优跨线传输率配置（3条线×8个缓冲区）：\n{np.round(final_TS, 3)}")
        return final_C, final_TS, final_best


# -------------------------- 4. 算法运行 --------------------------
if __name__ == "__main__":
    config = ProblemConfig()
    c_mohfa = C_MOHFA(config)
    final_C, final_TS, final_obj = c_mohfa.run(max_iter=100)
